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1. Introduction and General Status 

The main goals of the research under this grant consist of the development of 
mathematical tools and measurement of transport properties necessary for high fidelity modelling 
of crystal growth from the melt and solution, in particular for the Bridgman-Stockbarger growth 
of mercury cadmium telluride (MCT) and the solution growth of triglycine sulphate (TGS). Of 
the tasks desribed in detail in the original proposal, two remain to be worked on: 

- ! development of a spectral code for moving boundary problems, 

diffusivity measurements on concentrated and supersaturated TGS solutions. 

During this seventh half-year period, good progress has been made on these tasks. 

2. MCT Code development 

In the last six-monthly report we described our work on the solution of the coupled 
equations governing momentum (fluid flow) and heat transport together with a moving boundary. 
We reported that the chief problems met during this period have been associated with excessively 
long iteration times and that a Preconditioned Generalized Conjugate Residual (PGCR) scheme 
had been implemented for an irregular domain with a fixed non-planar boundary in an attempt to 
alleviate this problem and that we would apply it to the moving boundary if appropriate. 

We have since successfully implemented this scheme for the general moving boundary 
problem. Our method and the results of a careful comparison of our Bridgman code with other 
work is described in the attached paper which will be submitted for publication in the 
International Journal of Numerical Methods for Heat and Fluid Flow. In an attempt to 
demonstrate the portability of our code from the CRAY to workstations we have tested the 
performance of our code on a variety of machines. The results are also described in the paper. 
Work will continue on further development of the code to include a non-dilute dopant and to 
investigate the feasibility of incorporating radiation between the furnace and ampoule using a 
Finite Element approach. 

3. Diffusivity Measurements 

In the last report we had stated that the development of a new interferometric method to 
measure the diffusivities of supersaturated solutions had been completed. The method was 
described in that report as well as in earlier ones. A rectangular optical cell is initially filled with 
a solution of a certain concentration C. A solution of higher concentration C + AC is then 
injected at the bottom of the cell by means of a syringe. The resulting system is convectively 
stable, in one-dimensional terms, with the heavier solution below the lighter one and mixing 
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occurs by diffusion only. A Zygo Mark n Mach-Zehnder interferometer interfaced with a 
personal computer is used to follow the evolution of the concentration profile in the cell. At 
regular intervals the interferometric intensity profiles produced are stored and the advanced 
fringe analysis software ZAPPC is used to convert these into refractive index profiles which are 
proportional to the concentration. A numerical integration of these profiles yields the diffusivity. 

Several measurements were done with undersaturated and supersaturated solutions of 
sodium chloride and some of these results were shown in our last report. Our method compares 
well with earlier measurements of diffusivity [1,2] and has a reproducability of better than 1% 
for undersaturated solutions. However, in the supersaturated region the only other attempt to 
measure diffusivities was done by Myerson and co-workers and our values are significantly 
lower than the value reported by them for sodium chloride. Their method, as are most other 
methods, is based on the assumption of constant diffusivity in the concentration range of the 
solution pair used in each measurement. In the supersaturated region the diffusivity is a strong 
function of the concentration and this assumption is no longer valid. Our method does not make 
such an assumption and this is probably the reason for the discrepancy between our values and 
theirs. 

Diffusivity measurements were also made with TGS solutions at 25°C and these are 
shown in fig. 1. The crystalline TGS used to prepare the solutions were provided by Dr. Roger 
Kroes of NASA Marshall Space Flight Center, He also provided the original data from his 
diffusivity and refractive index measurements for undersaturated TGS solutions that were 
presented in graph form in his publications [3,4]. Solubility data for TGS solutions were also 
obtained from the literature [3,5,6]. The refractive index of supersaturated TGS solutions was 
needed for the diffusivity calculations and this was measured with a temperature controlled 
Milton Roy Abbe-3L refractometer. 

We have shown that our method can be successfully used to measure the diffusivities of 
supersaturated TGS solutions at ambient temperatures. We are currently performing 
measurements at concentrations between those shon in fig. 1, in order to obtain a reliable curve 
for the concentration dependence of TGS diffusivities at 25° C. In the near future we will be 
designing a temperature controlled apparatus to measure diffusivities at other temperatures as 
well. 
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4. Presentations and Publications 

From the work carried out under this grant the following papers have been published, 
accepted for publication or are in preparation for submission for publication: 

1. A. Nadarajah, F. Rosenberger and J. I, D. Alexander, Modelling the Solution Growth of 
Triglycine Sulfate in Low Gravity, J. Crystal Growth 104 (1990) 218-232. 

2. F. Rosenberger, J. I. D. Alexander, A. Nadarajah and J. Ouazzani, Influence of Residual 
Gravity on Crystal Growth Processes, Microgravity Sci. Technol. 3 (1990) 162-164. 

3. J. P. Pulicani and J. Ouazzani, A Fourier-Chebyshev Pseudo-Spectral Method for Solving 
Steady 3-D Navier-Stokes and Heat Equations in Cylindrical Cavities, Computers and 
Fluids 20 (1991) 93. 

4. J. P. Pulicani, S. Krukowski, J. I. D. Alexander, J. Ouazzani and F. Rosenberger, 
Convection in an Asymmetrically Heated Cylinder, Int. J. Heat Mass Transfer 35 (192) 
2119. 

5. F. Rosenberger, J. I. D. Alexander and W.-Q. Jin, Gravimetric Capillary Method for 
Kinematic Viscosity Measurements, Rev. Sci, Instr. 63 (192) 269. 

6. A. Nadarajah, F. Rosenberger and T, Nyce, Interferometric Technique for Diffusivity 
Measurements in (Supersaturated) Solutions, J. Phys. Chem (submitted). 

7. F. Rosenberger, Boundary Layers in Crystal Growth, Facts and Fancy, in Lectures on 
Crystal Growth, ed. by H. Komatsu (in print). 

8. F. Rosenberger, Short-duration Low-gravity Experiments - Time Scales, Challenges and 
Results, Microgravity Sci. Applic. (submitted), 

9. Y. Zhang, J.I.D. Alexander and J. Ouazzani, A Chebishev Collocation Method For 
Moving Boundaries, Heat Transfer and Convection During Directional Solidification, 
Intemat. J. Numerical Methods Heat Fluid Flow (in preparation). 

In addition to the above publications, the results of our work have been presented at the 
following conferences and institutions: 

1. J.I.D. Alexander, Modelling the Solution Growth of TGS Crystals in Low Gravity, 
Committee on Space Research (COSPAR) Plenary Meeting, The Hague, Netherlands, 
June 26 - July 6, 1990. 

A. Nadarajah, Modelling the Solution Growth of TGS Crystals in Low Gravity, Eighth 
American Conference on Crystal Growth, Vail, Colorado, July 15-21, 1990. 


2 . 



6 


3. J.I.D. Alexander, Commercial Numerical Codes: To Use or Not to Use, Is This The 
Question?, Microgravity Fluids Workshop, Westlake Holiday Inn, Cleveland Ohio, 
August 7-9, 1990. 

4. F. Rosenberger, Fluid Transport in Materials Processin, Microgravity Fluids Workshop, 
Westlake Holiday Inn, Cleveland Ohio, August 7-9, 1990. 

5. F. Rosenberger, Influence of Residual Gravity on Crystal Growth Processes, First 
International Microgravity Congress, Bremen, September 1990 (invited). 

6. J.I.D. Alexander, Residual Acceleration Effects on Low Gravity Experiments, Institute de 
Mechaniques des Fluides de Marseilles, University de Aix-Marseille HI, Marseille, 
France, January 1991, (3-lecture series, invited). 

7. J.I.D. Alexander, An Analysis of the Low Gravity Sensitivity of the Bridgman-Stockbarger 
Technique, Department of Mechanical Engineering at Clarkson University, April 1991 
(invited). 

8. A. Nadarajah, Measuring Diffusion Coefficients of Concentrated Solutions, Fifth Annual 
Alabama Materials Research Conference, Birmingham September 1991. 

9. A. Nadarajah, Modelling Crystal Growth Under Low Gravity, Annual Technical Meeting 
of the Society of Engineering Science, Gainesville, November 1991. 

10. J.I.D. Alexander, Vibrational Convection and Transport Under Low Gravity Conditions, 
Society of Engineering Science 28th Annual Technical Meeting, Gainesville, Florida, 
November 6-7 , 1991. 

11. F. Rosenberger, Theoretical Review of Crystal Growth in Space - Motivation and Results, 
International Symposium on High Tech Materials, Nagoya, Japan, November 6-9, 1991 
(plenary lecture, invited). 

12. F. Rosenberger, Computer Simulation in Materials Science, Mitsubishi Frontiers 
Research Institute, Tokyo, Japan, November 8, 1991 (invited). 

13. F. Rosenberger, Importance of Materials Research in Space Laboratories for Industrial 
Development, International Symposium for Promoting Applications and Capabilities of 
the Space Environment, Tokyo, Japan, November 14-15, 1991 (plenary lecture, invited). 

14. F. Rosenberger, What Can One Learn from 10 Second Low-Gravity Experiments?, In 
Space 1991, Tokyo, Japan, November 14-15, 1991 (plenary lecture invited). 

15. P. Larroude, J. Ouazzani and J.I.D. Alexander, Flow Transitions in a 2D Directional 
Solidification Model, 6th Materials Science Symposium, European Space Agency, 
Brussels, Belgium, 1992 (poster). 

16. F. Rosenberger, Microgravity Materials Processing and Fluid Transport, AIAA Course 
on Low-Gravity Fluid Dynamics, AIAA Meeting, Reno, NY, January 10-12, 1992 (3- 
lecture series, invited). 
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17. J.I.D Alexander, Numerical Simulation ofLow-g Fluid Transport, AIAA Course on Low- 
Gravity Fluid Mechanics, Reno, NV, January 10-12, 1992 (invited). 

18. F. Rosenberger, Time Scales in Transport Processes and Challenges for Short-Duration 
Low-Gravity Experiments, Falltower Days Bremen, Bremen, Germany, June 1-3, 1992 
(invited). 

19. J.I.D. Alexander, Modelling or Muddling? Analysis of Buoyancy Effects on Transport 
under Low Gravity Conditions, World Space Congress, Washington, DC, August 28 - 
September 5, 1992 (invited lecture). 



A CHEBYSHEV COLLOCATION METHOD FOR MOVING 
BOUNDARIES, HEAT TRANSFER, AND CONVECTION 
DURING DIRECTIONAL SOLIDIFICATION 
Yiqiang Zhang, J. Iwan D. Alexander and Jalil Ouazzani* 

Center for Micogravity and Materials Research, University of Alabama in Huntsville 

Abstract 

Free and Moving Boundary problems require the simultaneous solution of unknown field variables 
and the boundaries of the domains on which these variables are defined. There are many 
technologically important processes that lead to moving boundary problems associated with fluid 
surfaces and solid-fluid boundaries. These include crystal growth, metal alloy and glass 
solidification, melting and flame propagation. The directional solidification of semi-conductor 
crystals by the Bridgman-Stockbarger method 1 ' 2 is a typical example of a such a complex process. 
A numerical model of this growth method must solve the appropriate heat mass and momentum 
transfer equations and detemine the location of the melt-solid interface. In this work, a Chebyshev 
pseudospectral collocation method is adapted to the problem of directional solidification. 
Implementation method involves a solution algorithm that combines domain decomposition, a 
finite-difference preconditioned conjugate minimum residual method and a Picard type iterative 
scheme. 


* Presently at the Institute de Mecanique des Fluides de Marseille, 1 rue Honnorat, Marseille, France. 



2 


1. INTRODUCTION 

Moving and free boundary problems arc problems that require as part of the solution the 
determination of some or all the boundaries of the domain under consideration. Included in this 
class of problems are situations that involve fluid surfaces, or solid-fluid interfaces. Freezing and 
melting, crystal growth, flame propagation, liquid surface configurations, are examples of such 
processes that are important in a variety of areas with technological applications. Such problems 
generally pose a challenging problem to the numerical modeller. The Bridgman-Stockbarger 
directional solidification crystal growth technique is a typical example of such a complex problem. 
To adequately represent the physics of the problem, the solution method must be able to cope with 
the following: The unknown location of the crystal-melt interface, high Rayleigh number 
buoyancy-driven flows, heat transfer by conduction (along ampoule walls and in the crystal), 
convective-diffusive heat transfer in the melt and radiative and convective heat transfer between the 
furnace and the ampoule. Even for pure melts, due to differences in thermal conductivities between 
melt , crystal and ampoule, and the differences in thermal and momentum diffusivities in the melt, 
the problem has a variety of disparate length scales over which characteristic features must be 
accurately represented. 

In past work 3-11 , the Finite Element Method (FEM) has been successfully applied to the 
problem of computing melt and crystal temperature and concentration distributions, melt 
convection and the location of the crystal-melt interface. As an alternative to FEM we present a 
Chebyshev collocation(pseudospectral) method suitable for the solution of this class of problem. 
Spectral and pseudospectral methods 12 * 13 involve the representation of the solution as a truncated 
series of smooth functions of the independent variables. In contrast to FEM, for which the 
solution is approximated locally with expansions of local basis functions, spectral methods 
represent the solution as an expansion in global functions. In this sense they may be viewed as an 
extension of the separation of variables technique applied to complicated problems 14 . 

For problems that are characterized either by irregularly shaped domains, or even domains 
of unknown shape, it is, in general, neither efficient nor advantageous to try to find special sets of 
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spectral functions that are tuned to the particular geometry in consideration (especially in the case of 
solidification, where the melt-crystal geometry is not known a priori). Two alternative methods 
are mapping and patching 14 . Mapping allows an irregular region to be mapped into a regular one 
(which facilitates the use of known spectral functions, such as Chebyshev polynomials). For 
directional solidification systems (see Fig. 1) the melt-crystal boundary and, thus, the melt and 
crystal geometries, are unknown. Nevertheless, by specifying the melt-crystal boundary as some 
unknown single- valued function, the melt and crystal geometries can be mapped into simple ones 
by a smooth transformation. This mapping facilitates the use of Chebyshev polynomials to 
approximate the dependent variables in these new domains. 

As can be seen from Figs. 1 and 2, heat transfer to and in the ampoule wall must also be 
considered. To do this we employ patching by subdividing the system into four domains (crystal, 
melt and two ampoule domains), and transform these domains to domains with simple shapes. We 
then solve the resulting problems in each domain and solve the full problem in the complicated 
domain by applying suitable continuity conditions across any boundaries (real or artificial) between 
the domains. 

The formulation of the problem is outlined in section 2. The solution method is described in 
section 3. Our results are presented in section 4 and discussed in section 5. 

2. FORMULATION 

The vertical Bridgman-Stockbarger system is depicted in Fig. 2. A cylindrical ampoule 
with inner and outer diameters of 2Ro and 2 (Rq+R w ) contains melt and crystal. To grow the crystal 
the ampoule must be translated relative to a prescribed external temperature gradient. The objective 
of this model is to describe a steady growth process that, in reality, can be achieved between initial 
and terminal transients in sufficiently long ampoules. Toward this end a pseudo-steady state 
model 2 is employed that neglects the ends of the ampoule. The remainder of the ampoule is 
assumed to occupy a cylindrical computational region of length L. Ampoule translation is then 
accounted for by supplying a melt to the top of the computational space at a uniform velocity, and 
withdrawing crystal from the bottom at the same velocity. It is thus assumed that there is no lag 
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between the translation rate and the crystal’s growth velocity. Transport of heat from the furnace to 
the ampoule is modelled using a prescribed furnace temperature profile. The heat transfer from the 
furnace to the outer ampoule wall is governed by a heat transfer coefficient Bi(z). This is discussed 
in more detail later. The top and bottom of the ampoule are respectively assigned temperatures of 
Th and Tc (Th>Tc). 

The variables are cast in dimensionless form by using Ro, oljJRq, RqOCl, ocl/Rq 2 and 
Th - Tc, where ocl is the melt’s thermal diffusivity, to scale length, velocity, stream function, 
vorticity and temperature, respectively. That is, 

x = (r,z)=( r,z)/Ro, u = u R</cxl, \| f =\fl Rc«l , <*> = coRotac, T= J' - JS- • (1) 

Th-Tc 

Here r and z represent the radial and axial coordinates, \j/ is the stream function, to is the vorticity 
and u = (u,w) represents velocity with radial and axial components u and w, respectively. A tilde 
denotes a dimensional quantity. Melt, crystal and ampoule temperatures will be distinguished by 
the suffixes L (melt), S (crystal)and W (ampoule) when necessary. The location of the crystal melt 
boundary is given by z=h(r,t) and must be determined. The melt is assumed to be incompressible 
and the stream function and vorticity are defined by the velocity components (u,w) as 


u 


_i<ty 
r dz ’ 


-1 

w= T^ m 


du 

dz 


dw 
dr • 


The governing equations then take the following form 
In the melt, 0 <r < 1, 0< z < h(r,z) 


( 2 ) 


u §£ +w 3“_m=p r |^“+l?2. + &)_p r m. 

dr dz r \ dr 2 r dr 3z 2 ) i 2 


rco 


dr 2 r 9r 


a ^)- pr s +PAir ^ • 

dz 2 1 i 2 dr 

(3) 

8z 2 ’ 

(4) 


and 
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3T 9T A - . 3^19 3 2 

U ¥ +W ^ =AT ■ (5) 

where Pr = v/ocl is the Prandtl number, Gr = (3(Th- Tc)gRo 3 /v 2 is the Grashof number, v is the 

melt’s viscosity and (3 is the melt’s thermal expansion coefficient 

In the crystal, 0 < r < 1, h(r) « z < A , 

<x'Pe^=AT, (6) 

oz 

and in the ampoule wall, 1 < r < r w , Ck z < A, 

a"p<||=AT, (7) 

where a', and a" are, respectively, the ratios of the melt’s thermal diffusivity with the crystal and 
ampoule thermal diffusivities, and Pe = VoRo/ccl is the Pedet number and Vo is the ampoule 
translation rate. 

For the temperature the boundary conditions are: 


At the melt-crystal interface z= h(r,t) 


Tl=Ts = Tm, 


k' VT lii - VT sn=StPe a'ne z , 


where Tm represents the dimensionless melting temperature, k' is the ratio of melt and crystal 
conductivities, St = AH/(C ps AT) is the Stefan number. The vector n is the unit normal to the 
crystal-melt surface and points into the melt. At the outer ampoule wall, r= r w 


rfT 

|UBi(zXT-Tp(z)). 


The temperatures at z=0 and z=A are constant, i.e. 


T(r,0) = 1, T(r,A) = 0, 
and the heat flux is continuous across the inner ampoule wall 


( 11 ) 
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aT(i,z)\ _ l j dm,z) \ 
dr )l \ dr )w 



In (10) Bi(z) is a heat transfer coeffcient and Tf(z) is the furnace temperature profile. The 
coefficients k* and k** represent the ratio of the wall conductivity with that of the melt and 
ampoule, respectively. 

For the stream function the boundary conditions are 

y (0,z) = 0, Y(l,z) = Pe, y(0,z) = - i r 2 Pe, \j/(h(r),z) = - ^ r 2 Pe, (13) 

and the vorticity is zero at r=0. At the other melt boundaries the boundary conditions for the 
vorticity are enforced (iteratively) using previously computed values of the velocity field (the 
scheme is explained in section 3.3. The velocity boundary conditions are 

u(0, z) = u(l,z)=u(r,h(r)) =0; w(0,z) = w(l,z)= w(r,h(r))=Pe. (14) 

Note that, at the melt-crystal boundary there are two boundary conditions for the 
temperature. In the following section we describe an iterative scheme which distinguishes one of 
the temperature boundary conditions and uses it to compute the interface shape iteratively. 

3. SOLUTION METHOD 

The solution method is based on a Picard 15 type iteration which consists essentially of four 

steps: 

1. The initial shape of the crystal-melt interface is specified and an independent variable 
transformation is applied to the governing equations and boundary conditions in the melt, crystal 
an ampoule regions. This specifies the computational domains. 

2. The coupled momentum, heat, mass and species equations are then solved using three of the 
four boundary conditions on the moving boundary. 
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3. The remaining boundary condition (or distinguished condition 2 ), in this case equation ( 8 ), is 
used to compute corrected boundary locations. 

4. Steps 2 and 3 are repeated until the distinguished boundary condition is satisfied. 

The solution method is implemented using domain decomposition and a preconditioned 
generalized conjugate residual method 13 ’ 16 
3.1 Domain Decomposition 

The physical region is split into four computational domains, Qi, i=l,...,4. The domains 
correspond to the melt (Si), the crystal (S 3 ), and the portions of the ampoule wall adjacent to the 
melt (S 2 ) and the crystal (S 4 ). The irregularly shaped domains are mapped onto rectangular 
regions by 


l=r ’ n= hlr £ ^ Qi ' 

(15) 

^ =r,T1= ha7 ,E2 '^ 

(16) 


(17) 

^ =r , T| = 2 - h(1) _ A , - 4-> n 4. 

(18) 


3.2 Spatial discretization 

The dependent variables, O are approximated by Chebyshev polynomials 12 ’ 13 , i.e 

N.M 

<I>=Onm (Xi,Yj)= X a ijTij(Xi,Yj ), (19) 

i = 0 

j - 0 

where Ty = TjTj , and the Tk are Chebyshev polynomials of order k. The points (X{ , Yj ) are 

related to the coordinates \ and T| by 

^ = aX + b, T] = c Y + d, ( 20 ) 

where a and b are determined by the transformation of each domain, Qj, to [-l,l]x[— 1 , 1 ]. The 
discrete points (Xj,Yj), i=0, N, j=0,M,are the Gauss-Lobatto collocation points 13 . That is, 
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X j = cos x 


1 

N 


,i = 0, 1„..,N 


( 21 ) 


X ; = cos n 


l 

N 


,i = 0, 


( 22 ) 


The spatial derivatives are given by 


30 _ i 80 80 _ i 80 3 2 0 _ i 8 z O 8 2 0_ i 8 2 o 3 2 0_ j 8 2 o 
9 ^ a 9 X’ 9 r) cay’ ac 9 x 3 y’ 9^2 a 2 0x 2’ 9^2 c 2 ay 2 

where the derivatives with respect to X and Y have the forms 

£ D ^ (x ^ Yj)= | D .P 0pj , 

3X 5X p = 0 p = 0 

»(X.Y).i»(X t Y J >_ | ^ (XiiYq)= £ 

9Y 8Y q = 0 q = 0 

£ D&0(X p ,Y j )» £ D^pi, 

8a 3a p=o p=o 


(23) 


(24) 


(25) 


(26) 


£ Dy q y <D (Xj, Yq ) = £ D^i,, 

dx dx 


q = 0 


q = 0 


tftolY Y i N ’ M N ’ M 

= X I^x P D , y q ^(Xp,Y q )= X Dj^O 

9XY p = 0 


p=0 

q = 0 q = 0 


T* pq« 


(27) 


( 28 ) 


where the expressions for D x , D y , Dxx, D yy and Dxy are given explicitly by Ouazzani 17 . 


3.3 Pseudo-unsteady discretization 
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The governing equations now have the form 

(A® - (i) = S« i = 1,6 (29) 

where the A®, A(0, S(>) and <j>ti) are given in Appendix I. To solve these equations using a 
pseoudo- unsteady iterative scheme (28) is rewritten as 

+ A® »> (i) - S® + A (i fy (i) , i = 1,6 (30) 

dx 

The left-hand side of (29) is written in discrete form as 

(^ + A)(j>={A + aW +1 - o<t> n , (31) 

where o(0 = 1/Ax® for i=l,2,4,5,6 and is zero for i=3. the index in parentheses has been omitted 
for clarity and the superscript denotes the pseudo-time or iterative step number. Note that the time 
step size, At®, is generally different for each of the equations. 

The problem now has the form 

H sp (|) n+1 + oI<j> n+1 = F(<|),h) n , (32) 

where 

H sp = A -A, and F= S + G<J> n , (33) 

and H sp j s obtained from the expressions in Appendix B using the Chebyshev derivatives (26)- 
(28) and equation (23). A superscript n denotes a quantity evaluated at the nth iterative step (note 
that the indices in parethenses have been omitted for clarity) . 

3.4 Vorticity boundary condition 

To solve the vorticity-stream function equations we adopted the following procedure which 
is simply an extension, for Chebyshev approximations, of an approach described by Peyret 12 . 
The velocity field is calculated from the stream function obtained from the previous iteration. The 
vorticity at the boundary which corresponds to this velocity field is then found from 
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~n+l_$U dwf 

“\3z 3r I’ 

and the value of the vorticity to be applied at the boundary, 0> n+1 , is given by 

©n+l =^ n+1 + (1 — y)0) n . 

Here y(0<Y<l)isa relaxation parameter. 

3.5 Preconditioned Generalized Conjugate Residual Method 
The operator H sp is represented by a full matrix of order (N+l) 2 x (M+l) 2 and is not 
symmetric. In order to solve the system of equations and boundary conditions represented by 
(28)-(32) and (A.9)-(A.18), each of the spectral operators H sp for each of the domains Qi , i=l,4 
and the conditions on shared domain boundaries Qi n i*j, i,j=l,4 are combined and 
approximated by a single finite difference operator Hfd. This is defined over the entire domain 
=^X2j. The following iterative procedure which consisits of inner and outer loops is then 

adopted: 

Outer loop: First an initial interface shape h° is assumed 
Inner loop: The residual R is then initialized by 

R° = H S ;®-F, 

where O represents the <j)®. Then we solve 

Hfd®° = R°» 

where H* = H + cl. Then we set 


(36) 

(37) 


( 34 ) 


(35) 


P° = ©°, 


(38) 


and calculate 



( 39 ) 


(R m ,H sp P m ) 
(H sp P In »H S pP m ) ' 


The variables O are then updated from 


O m+ i = O m + a^P" 1 , 


(40) 


and the problem 


Hf d 0 m+1 = R m+1 , 


(41) 


is solved for 0. P is then updated using 


p"* 1 = © m + Z Bf fl P j , 
j— o J 

where 

pnH-i _ (Hsp 0 ^ 

3 (H sp pi,H sp pi) ’ 


(42) 


(43) 


The procedure is continued until | R | < e. 

The preconditioned problem is given by equations (37) and (41). The finite difference 
operator H*fd is approximated by incomplete LU decomposition. The solution for 0 is obtained 
by forward and backward substitution. The subsequent approximations to O = (Ts,Tl,£0,\|/ ) are 
then obtained from (40). At this point we note that while we used a nine-diagonal matrix for the 
second-order central finite difference operator for the solution of the temperature field, a seven 
diagonal operator was used for the solution of the stream-function and vorticity as it appeared to 
lead to more rapid convergence. This means that the cross-derivative terms were evaluated at the 
previous time step and were included in F on the right-hand side of (32). 

3.6 Interface Shape Update 



12 


This iterative procedure is repeated until the convergence criterion is satsified. The first of 
equations (A. 17) is used as a distinguished boundary condition. If it is not satisfied, another outer 
loop iteration is performed and the interface shape is relocated using either Newton’s method 



where 0i is the difference between the temperature at the ith interfacial site and the melting 
temperature T m ; or from a searching method 

h? +1 = h? + aCr r T M ). (45) 

Here a is found by numerical experiment. We found that by using the Newton method for the first 
few iterations and then the searching method for subsequent iterations, we achieved better success 
than with the Newton method alone. 

4, RESULTS 

We carried out several tests of the method. The results are shown in Table 1 and in Figs. 3 
and 4. The parameters used are given in Appendix C and correspond to the thermophysical 
properties of Gallium-doped Germanium. For the cases examined our results are in good 
agreement with the FEM calculations of Adomato and Brown. 2 For example, the difference 
between our results and theirs for the maximum stream function computed for Fig. 3c is less than 
2 %. 

Figure 2 shows results for a furnace with a constant temperature gradient and Bi=7.143. 

That is, 

Tf(z) = 1 - zA' 1 . (46) 

The isotherms are practically flat except at the crystal-melt boundary where the mismatch in thermal 
conductivity results in a convex interface. The flow depicted by the streamlines in Fig. 3b-d is a 
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combination of the ampoule translation (which, if buoyant convection were absent, would appear 
as a set of vertical streamlines parallel to the ampoule wall) and buoyant flow caused by radial 
gradients in temperature. This results in an upward flow of hot melt near the ampoule wall and a 
downflow near the ampoule centerline. Note the increase in flow intensity as the Grashof number 
is increased. 

Figure 4 shows results for different Grashof numbers for a non-uniform furnace 
temperature 
profile 

Tf(z) =0.5[1+ tanh (6-12Z/T 1 )] (47) 

together with a position dependent heat transfer function given by 

Bi(z) =0.2{2[1+ tanh (5-2z)] + 1+ tanh (2z-15)}. (48) 

Radial temperature gradients arise for two reasons in this problem: The mismatch in thermal 
conductivities at the ampoule-melt-crystal junction and the change in heat transfer at the quasi- 
adiabatic zones. These zones are created by the furnace temperature profile and conditions (47) and 
(48). This heating configuration produces two counter rotating cells. The upper cell increases in 
spatial extent as the Grashof number is increased. 

Table 1 shows the CPU times, number of iterations taken to converge and compiler options 
for the case shown in Fig. 4b for a CRAY/XMP, an iPSC parallel processor and an Ardent Titan 
computer. 

5. DISCUSSION 

Chebyshev spectral methods that have been shown to achieve superior accuracy for a wide 
range of fluid flow problems defined in regular geometries can be applied to problems involving 
unknown free and moving irregular boundaries through a combination of mapping and domain 
decomposition. For the directional solidification described here, this was achieved without 
incurring excessive CPU times and has been implemented on several different machines to 
illustrate the magnitude of the CPU times involved for a typical calculation. Whether there is 
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ultimately any advantage in using such spectral methods over finite elements will depend on the 
specific application. It will most likely depend on the accuracy required and on whether the ability 
of the Chebyshev collocation method to achieve better accuracy for a given number of collocation 
points (which is recognized for a variety of flows in regular geometries) is retained or degraded 
when using domain decomposition. 
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Appendix A 
Transformed Equations 

After the equations and boundary conditions (2) - (14) have been transformed according to 
(15)-(18) we obtain the following equations. 

For 0 < tj < 1 

0 < £ < 1 


where 



1 BT 9v|/ _ 3T dy y _ A * T 

^ an ar 

(A.1) 

l ByBcO 

on 


(A.2; 



(A.3) 



(A.4) 


and 
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A = J^+P-dhf B = ~ _jjl dh c ^ [2fidh]2_Id3i — l_dh 
h 2 lh dr/ ’ h dr ’ ^ dr/ Eh dr 


di 2 Ejh 


For 1 < tj < 2 
0 < ^ < 1 

where 

and 


A**T-Pea'i^- = 0, 

h^n 


B 2 3 2 
A** = — — + A*— — + B* 


a 2 + iA +c *A 


3^ 2 an 2 a$an 4 34 an 

A* l_ + (j!^dkf,B»=-=2-(n-2) 

rh-A^ \h - A dr; h-A dr 


(h-A) 2 

C* =(i 


W 1 dh\ 2 

1 d% 

1 dh 1 

\h — A dr 1 

h“ A di 2 

^(h-A) dr 


(A.5) 

(A.6) 

(A.7) 

(A. 8) 


In the ampoule wall, 1 < E, < r w , 0 < H <2, where h(r) is taken to be a constant at each inner 
iteration, we have 


and 


a*r t ! a*r 
5E, 2 hW 


+ 


i 5 I 

43| 


- pea if=°' 


0 <T1 < 1, 


S *T , 1 3^ , 

a? (h-A]?9n 2 


% 34 


'+Peoc 


- 1 3T 
h- A^H 


= 0, 1 <rj < 2. 


The boundary conditions become 

— = 0, \j/ = 0, © =0 at % = 0, 
OS 


(A.9) 
(A. 10) 
(A. 11) 


T = 1, y = - l/2^Pe, arn = 0, 


(A. 12) 
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Hcs 

1 

II 

> 



(A. 13) 


(A. 14) 

Hr = Bi(n)(T - T f (tj)), % = r w , 0 < T| <2 , 

(A. 15) 

O 

ft 

li 

N> 

O 

A 

xTr* 

A 

I 4 

(A. 16) 

Finally, at the crystal melt interface the boundary conditions are 


T = T m , l k^ L + 1 ?? = StPea ' , 

h dn h-A <h i^dhj 2 

(A. 17) 

and 


y = -^Pe, atrj = 1,0<£< 1 . 

(A. 18) 

In (A. 17) we have used the fact that the melting temperature Tm is assumed to be constant along 
the crystal melt interface (i.e. 8T/8£ = 0). The vorticity boundary condition is given by equation 

(35) with 

£_1 du t\ dh dw 8w 
h 8r| h dr 8r| 8£ 

(A. 19) 

Appendix B 


The Ati), Ati) and F® referred to in section 3.3 are expressed in terms of the equations 

given in Appendix I as follows: 
4>(l) = T n+1 


.<» _i ^T n+1 ar OT" +1 ars 

3| 8q drj 8£ 

(B.l) 

A (1) = A* 

(B.2) 
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ptt) = oW 

(B.3: 

<j>(2) = ojn+1 

A (2^ 1 /^V" 3tO n+1 dy 11 9w n+1 ojn+1 3y n . 

§i3n 3| 35 3n 5 3n ’ 

(B.4) 


A<» =A*- Et 

S 2 

(B.5) 


FO>=P^J 3 ^‘^* 3T " +1 ) + o<2V 0 « 
\ % h dr an / 

(B.6; 

<j>(3) = \j/n+l 

a< 3 > =0 

(B.7) 


A (3) =a*-2JL, 

m 

(B.8) 


f ( 3 > = £co n+1 , 

(B.9) 

<|>( 4 ) = T w n+1 

(0<n<l), <t>< 5 > = Tw n+1 (1 <ri < 2), <j)(6) = T s n+1 , 



A< 4 > = A< 5 > =A< 6 ) =0 

(B.10) 


a ( 4 ) a 2 1 a 2 1 a n ,a a 

A w = —t+- l ~- ■■ ■ + 4-— - Pea 

a ^ 2 h 2 an 2 h ^n 

(B.n; 


c 5 ) s 2 1 a 2 1 a „ 1 a 

A (5) =^-+ — 1 „ ■+ *■“ +Pea - 5 — 

a ^ 2 (h-A^an 2 h-Aan 

(B. 12 ; 


A<6) = A** Pe a' l J* , 

han 

(B.13) 


F® = oa^® , i = 4,5,6 

(B.14) 
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Appendix C 

Physical constants, system dimensions and thermophysical properties of Gallium doped 
Germanium used in the calculations 


Property 

dimension 

Ge:Ga 

Growth velocity 

[cms- 1 ! 

4X10- 4 

Ampoule length (L) 

[cm] 


Constant gradient furnace (Fig. 2) 


7.0 

Heat pipe furnace (Fig. 3) 


7.62 

Outer ampoule radius (R w ) 

[cm] 


Constant gradient furnace 


0.7 

Heat pipe furnace 


0.952 

Inner ampoule radius (Rq) 



Constant gradient furnace 


0.5 

Heat pipe furnace 


0.762 

Kinematic viscosity (v) 

[cm 2 s' 1 ] 

h- * 

La 

/-N 

o 

'W' 

■ 

Thermal conductivity (ampoule) 

[WK-W 1 ] 


Constant gradient furnace 


3.27 

Heat pipe furnace 


0.26 

Thermal conductivity (crystal) 

[WK-W 1 ] 

0.17 

Thermal conductivity (melt) 

[WK-fcnr 1 ] 

0.39 

Density (crystal) 

[gcm-3] 

5.5 

Density (melt) 

[g cm* 3 ] 

5.5 

Heat of solidification (AH) 

Os- 1 ! _ 

460 

Specific heat (melt) 

JK-lg- 1 .... 

0.39 

Specific heat (crystal) 

JK-lg-l 

0.39 

Thermal expansion coefficient 

[K-l] 

5 (10)- 4 
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Figure Captions 

Fig. 1 Typical Bridgman-Stockbarger set-up 

Fig. 2 a) The model Bridgman-Stockbarger system and b) the computational domains 
Fig.3 Results for results for a furnace with a constant temperature gradient, Bi=7. 143 and a Pr = 
0.07 melt, a) Gr = 5206, b) Gr = 52,060 c) Gr = 520,600 
Fig. 4 Results for a non-uniform furnace temperature profile (47) and position dependent heat 
transfer coefficient (48) for Pr =0.007 and a) Gr = 7,140, b) Gr = 14280 c) Gr=7 1,400 
d) Gr= 142,800 . 
















FIG. 4 


machine 

compiler 

options 

cpu 

time 

(sec) 

number 

of 

outer 

iterations 

number of inner iterations 
for each outer loop 






cray 

XMP/24 

vector, double 
precision 

1157 

10 

1501 

340 

247 

190 

174 

150 

111 

87 

103 

94 

y 

vector, double 
precision and 
maximum 
optimization 

8076 

10 

613 

1342 

268 

275 

236 

174 

156 

151 

29 

65 

Ardent 
Titan II 
2 CPU’s 

vector, parallel, 
double 
precision and 
maximum 
optimization 

11015 

10 

613 

1342 

268 

275 

236 

174 

156 

151 

29 

65 


Table 1 : Comparison of CPU times for different machines for the Gr = case using our method. 














































